The textbook for the Data Science course series is freely available online.

Learning Objectives

How linear regression was originally developed by Galton
What confounding is and how to detect it
How to examine the relationships between variables by implementing linear regression in R

Course Overview

There are three major sections in this course: introduction to linear regression, linear models, and confounding.

Introduction to Linear Regression

In this section, you’ll learn the basics of linear regression through this course’s motivating example, the data-driven approach used to construct baseball teams. You’ll also learn about correlation, the correlation coefficient, stratification, and the variance explained.

Linear Models

In this section, you’ll learn about linear models. You’ll learn about least squares estimates, multivariate regression, and several useful features of R, such as tibbles, lm, do, and broom. You’ll learn how to apply regression to baseball to build a better offensive metric.

Confounding

In the final section of the course, you’ll learn about confounding and several reasons that correlation is not the same as causation, such as spurious correlation, outliers, reversing cause and effect, and confounders. You’ll also learn about Simpson’s Paradox.

Section 1 - Introduction to Regression Overview

In the Introduction to Regression section, you will learn the basics of linear regression.

After completing this section, you will be able to:

Understand how Galton developed linear regression.
Calculate and interpret the sample correlation.
Stratify a dataset when appropriate.
Understand what a bivariate normal distribution is.
Explain what the term variance explained means.
Interpret the two regression lines.

This section has three parts: Baseball as a Motivating Example, Correlation, and Stratification and Variance Explained.

Motivating Example: Moneyball

The corresponding section of the textbook is the case study on Moneyball

Key points

Bill James was the originator of the sabermetrics, the approach of using data to predict what outcomes best predicted if a team would win.

Baseball basics

The corresponding section of the textbook is the section on baseball basics

Key points

The goal of a baseball game is to score more runs (points) than the other team.
Each team has 9 batters who have an opportunity to hit a ball with a bat in a predetermined order.
Each time a batter has an opportunity to bat, we call it a plate appearance (PA).
The PA ends with a binary outcome: the batter either makes an out (failure) and returns to the bench or the batter doesn’t (success) and can run around the bases, and potentially score a run (reach all 4 bases).
We are simplifying a bit, but there are five ways a batter can succeed (not make an out):

Bases on balls (BB): the pitcher fails to throw the ball through a predefined area considered to be hittable (the strike zone), so the batter is permitted to go to first base.
Single: the batter hits the ball and gets to first base.
Double (2B): the batter hits the ball and gets to second base.
Triple (3B): the batter hits the ball and gets to third base.
Home Run (HR): the batter hits the ball and goes all the way home and scores a run.

Historically, the batting average has been considered the most important offensive statistic. To define this average, we define a hit (H) and an at bat (AB). Singles, doubles, triples and home runs are hits. The fifth way to be successful, a walk (BB), is not a hit. An AB is the number of times you either get a hit or make an out; BBs are excluded. The batting average is simply H/AB and is considered the main measure of a success rate.

Bases on Balls or Stolen Bases?

The corresponding section of the textbook is the base on balls or stolen bases textbook section

Key points

The visualization of choice when exploring the relationship between two variables like home runs and runs is a scatterplot.

Code: Scatterplot of the relationship between HRs and runs

if(!require(Lahman)) install.packages("Lahman")

## Loading required package: Lahman

if(!require(tidyverse)) install.packages("tidyverse")

## Loading required package: tidyverse

## ── Attaching packages ─────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

if(!require(dslabs)) install.packages("dslabs")

## Loading required package: dslabs

library(Lahman)
library(tidyverse)
library(dslabs)
ds_theme_set()

Teams %>% filter(yearID %in% 1961:2001) %>%
    mutate(HR_per_game = HR / G, R_per_game = R / G) %>%
    ggplot(aes(HR_per_game, R_per_game)) + 
    geom_point(alpha = 0.5)

Code: Scatterplot of the relationship between stolen bases and runs

Teams %>% filter(yearID %in% 1961:2001) %>%
    mutate(SB_per_game = SB / G, R_per_game = R / G) %>%
    ggplot(aes(SB_per_game, R_per_game)) + 
    geom_point(alpha = 0.5)

Code: Scatterplot of the relationship between bases on balls and runs

Teams %>% filter(yearID %in% 1961:2001) %>%
    mutate(BB_per_game = BB / G, R_per_game = R / G) %>%
    ggplot(aes(BB_per_game, R_per_game)) + 
    geom_point(alpha = 0.5)

Assessment - Baseball as a Motivating Example

What is the application of statistics and data science to baseball called?

[ ] A. Moneyball
[X] B. Sabermetrics
[ ] C. The “Oakland A’s Approach”
[ ] D. There is no specific name for this; it’s just data science.

Which of the following outcomes is not included in the batting average?

[ ] A. A home run
[X] B. A base on balls
[ ] C. An out
[ ] D. A single

Why do we consider team statistics as well as individual player statistics?

[X] A. The success of any individual player also depends on the strength of their team.
[ ] B. Team statistics can be easier to calculate.
[ ] C. The ultimate goal of sabermetrics is to rank teams, not players.

You want to know whether teams with more at-bats per game have more runs per game. What R code below correctly makes a scatter plot for this relationship?

Teams %>% filter(yearID %in% 1961:2001 ) %>%
    mutate(AB_per_game = AB/G, R_per_game = R/G) %>%
    ggplot(aes(AB_per_game, R_per_game)) + 
    geom_point(alpha = 0.5)

[ ] A.

Teams %>% filter(yearID %in% 1961:2001 ) %>%
    ggplot(aes(AB, R)) + 
    geom_point(alpha = 0.5)

[X] B.

Teams %>% filter(yearID %in% 1961:2001 ) %>%
    mutate(AB_per_game = AB/G, R_per_game = R/G) %>%
    ggplot(aes(AB_per_game, R_per_game)) + 
    geom_point(alpha = 0.5)

[ ] C.

Teams %>% filter(yearID %in% 1961:2001 ) %>%
    mutate(AB_per_game = AB/G, R_per_game = R/G) %>%
    ggplot(aes(AB_per_game, R_per_game)) + 
    geom_line()

[ ] D.

Teams %>% filter(yearID %in% 1961:2001 ) %>%
    mutate(AB_per_game = AB/G, R_per_game = R/G) %>%
    ggplot(aes(R_per_game, AB_per_game)) + 
    geom_point()

What does the variable “SOA” stand for in the Teams table?

Hint: make sure to use the help file (?Teams).

[ ] A. sacrifice out
[ ] B. slides or attempts
[X] C. strikeouts by pitchers
[ ] D. accumulated singles

Load the Lahman library. Filter the Teams data frame to include years from 1961 to 2001. Make a scatterplot of runs per game versus at bats (AB) per game.

Teams %>% filter(yearID %in% 1961:2001) %>%
  mutate(AB_per_game = AB / G, R_per_game = R / G) %>%
  ggplot(aes(AB_per_game, R_per_game)) + 
  geom_point(alpha = 0.5)

Which of the following is true?

[ ] A. There is no clear relationship between runs and at bats per game.
[X] B. As the number of at bats per game increases, the number of runs per game tends to increase.
[ ] C. As the number of at bats per game increases, the number of runs per game tends to decrease.

Use the filtered Teams data frame from Question 6. Make a scatterplot of win rate (number of wins per game) versus number of fielding errors (E) per game.

Teams %>% filter(yearID %in% 1961:2001) %>%
  mutate(win_rate = W / G, E_per_game = E / G) %>%
  ggplot(aes(win_rate, E_per_game)) + 
  geom_point(alpha = 0.5)

Which of the following is true?

[ ] A. There is no relationship between win rate and errors per game.
[ ] B. As the number of errors per game increases, the win rate tends to increase.
[X] C. As the number of errors per game increases, the win rate tends to decrease.

Use the filtered Teams data frame from Question 6. Make a scatterplot of triples (X3B) per game versus doubles (X2B) per game.

Teams %>% filter(yearID %in% 1961:2001) %>%
  mutate(doubles_per_game = X2B / G, triples_per_game = X3B / G) %>%
  ggplot(aes(doubles_per_game, triples_per_game)) + 
  geom_point(alpha = 0.5)

Which of the following is true?

[X] A. There is no clear relationship between doubles per game and triples per game.
[ ] B. As the number of doubles per game increases, the number of triples per game tends to increase.
[ ] C. As the number of doubles per game increases, the number of triples per game tends to decrease.

Correlation

The corresponding textbook section is Case Study: is height hereditary?

Key points

Galton tried to predict sons’ heights based on fathers’ heights.
The mean and standard errors are insufficient for describing an important characteristic of the data: the trend that the taller the father, the taller the son.
The correlation coefficient is an informative summary of how two variables move together that can be used to predict one variable using the other.

Code

# create the dataset
if(!require(HistData)) install.packages("HistData")

## Loading required package: HistData

library(HistData)
data("GaltonFamilies")
set.seed(1983)
galton_heights <- GaltonFamilies %>%
  filter(gender == "male") %>%
  group_by(family) %>%
  sample_n(1) %>%
  ungroup() %>%
  select(father, childHeight) %>%
  rename(son = childHeight)

# means and standard deviations
galton_heights %>%
    summarize(mean(father), sd(father), mean(son), sd(son))

## # A tibble: 1 x 4
##   `mean(father)` `sd(father)` `mean(son)` `sd(son)`
##            <dbl>        <dbl>       <dbl>     <dbl>
## 1           69.1         2.55        69.2      2.71

# scatterplot of father and son heights
galton_heights %>%
    ggplot(aes(father, son)) +
    geom_point(alpha = 0.5)

Correlation Coefficient

The corresponding textbook section is the correlation coefficient

Key points

The correlation coefficient is defined for a list of pairs \((x_1, y_1), ..., (x_n, y_n)\) as the product of the standardized values: \((\frac{x_i - \mu_x}{\sigma_x})(\frac{y_i -\mu_y}{\sigma_y})\).
The correlation coefficient essentially conveys how two variables move together.
The correlation coefficient is always between -1 and 1.

Code

rho <- mean(scale(x)*scale(y))

galton_heights <- GaltonFamilies %>% filter(childNum == 1 & gender == "male") %>% select(father, childHeight) %>% rename(son = childHeight)
galton_heights %>% summarize(r = cor(father, son)) %>% pull(r)

## [1] 0.5007248

Sample Correlation is a Random Variable

The corresponding textbook section is Sample correlation is a random variable

Key points

The correlation that we compute and use as a summary is a random variable.
When interpreting correlations, it is important to remember that correlations derived from samples are estimates containing uncertainty.
Because the sample correlation is an average of independent draws, the central limit theorem applies.

Code

# compute sample correlation
R <- sample_n(galton_heights, 25, replace = TRUE) %>%
    summarize(r = cor(father, son))
R

##           r
## 1 0.4787613

# Monte Carlo simulation to show distribution of sample correlation
B <- 1000
N <- 25
R <- replicate(B, {
    sample_n(galton_heights, N, replace = TRUE) %>%
    summarize(r = cor(father, son)) %>%
    pull(r)
})
qplot(R, geom = "histogram", binwidth = 0.05, color = I("black"))

# expected value and standard error
mean(R)

## [1] 0.4970997

sd(R)

## [1] 0.1512451

# QQ-plot to evaluate whether N is large enough
data.frame(R) %>%
    ggplot(aes(sample = R)) +
    stat_qq() +
    geom_abline(intercept = mean(R), slope = sqrt((1-mean(R)^2)/(N-2)))

Assessment - Correlation

While studying heredity, Francis Galton developed what important statistical concept?

[ ] A. Standard deviation
[ ] B. Normal distribution
[X] C. Correlation
[ ] D. Probability

The correlation coefficient is a summary of what?

[X] A. The trend between two variables
[ ] B. The dispersion of a variable
[ ] C. The central tendency of a variable
[ ] D. The distribution of a variable

Below is a scatter plot showing the relationship between two variables, x and y.

Scatter plot relationship x and y

From this figure, the correlation between x and y appears to be about:

[X] A. -0.9
[ ] B. -0.2
[ ] C. 0.9
[ ] D. 2

Instead of running a Monte Carlo simulation with a sample size of 25 from our 179 father-son pairs, we now run our simulation with a sample size of 50.

Would you expect the mean of our sample correlation to increase, decrease, or stay approximately the same?

[ ] A. Increase
[ ] B. Decrease
[X] C. Stay approximately the same

Instead of running a Monte Carlo simulation with a sample size of 25 from our 179 father-son pairs, we now run our simulation with a sample size of 50.

Would you expect the standard deviation of our sample correlation to increase, decrease, or stay approximately the same?

[ ] A. Increase
[X] B. Decrease
[ ] C. Stay approximately the same

If X and Y are completely independent, what do you expect the value of the correlation coefficient to be?

[ ] A. -1
[ ] B. -0.5
[X] C. 0
[ ] D. 0.5
[ ] E. 1
[ ] F. Not enough information to answer the question

Load the Lahman library. Filter the Teams data frame to include years from 1961 to 2001.

What is the correlation coefficient between number of runs per game and number of at bats per game?

Teams_small <- Teams %>% filter(yearID %in% 1961:2001)
cor(Teams_small$R/Teams_small$G, Teams_small$AB/Teams_small$G)

## [1] 0.6580976

Use the filtered Teams data frame from Question 7.

What is the correlation coefficient between win rate (number of wins per game) and number of errors per game?

cor(Teams_small$W/Teams_small$G, Teams_small$E/Teams_small$G)

## [1] -0.3396947

Use the filtered Teams data frame from Question 7.

What is the correlation coefficient between doubles (X2B) per game and triples (X3B) per game?

cor(Teams_small$X2B/Teams_small$G, Teams_small$X3B/Teams_small$G)

## [1] -0.01157404

Anscombe’s Quartet/Stratification

There are three links to relevant sections of the textbook:

Key points

Correlation is not always a good summary of the relationship between two variables.
The general idea of conditional expectation is that we stratify a population into groups and compute summaries in each group.
A practical way to improve the estimates of the conditional expectations is to define strata of with similar values of x.
If there is perfect correlation, the regression line predicts an increase that is the same number of SDs for both variables. If there is 0 correlation, then we don’t use x at all for the prediction and simply predict the average \(\mu_y\). For values between 0 and 1, the prediction is somewhere in between. If the correlation is negative, we predict a reduction instead of an increase.

Code

# number of fathers with height 72 or 72.5 inches
sum(galton_heights$father == 72)

## [1] 8

sum(galton_heights$father == 72.5)

## [1] 1

# predicted height of a son with a 72 inch tall father
conditional_avg <- galton_heights %>%
    filter(round(father) == 72) %>%
    summarize(avg = mean(son)) %>%
    pull(avg)
conditional_avg

## [1] 71.83571

# stratify fathers' heights to make a boxplot of son heights
galton_heights %>% mutate(father_strata = factor(round(father))) %>%
    ggplot(aes(father_strata, son)) +
    geom_boxplot() +
    geom_point()

# center of each boxplot
galton_heights %>%
    mutate(father = round(father)) %>%
    group_by(father) %>%
    summarize(son_conditional_avg = mean(son)) %>%
    ggplot(aes(father, son_conditional_avg)) +
    geom_point()

## `summarise()` ungrouping output (override with `.groups` argument)

# calculate values to plot regression line on original data
mu_x <- mean(galton_heights$father)
mu_y <- mean(galton_heights$son)
s_x <- sd(galton_heights$father)
s_y <- sd(galton_heights$son)
r <- cor(galton_heights$father, galton_heights$son)
m <- r * s_y/s_x
b <- mu_y - m*mu_x

# add regression line to plot
galton_heights %>%
    ggplot(aes(father, son)) +
    geom_point(alpha = 0.5) +
    geom_abline(intercept = b, slope = m)

Bivariate Normal Distribution

There is a link to the relevant section of the textbook: Bivariate normal distribution (advanced)

Key points

When a pair of random variables are approximated by the bivariate normal distribution, scatterplots look like ovals. They can be thin (high correlation) or circle-shaped (no correlation).
When two variables follow a bivariate normal distribution, computing the regression line is equivalent to computing conditional expectations.
We can obtain a much more stable estimate of the conditional expectation by finding the regression line and using it to make predictions.

Code

galton_heights %>%
  mutate(z_father = round((father - mean(father)) / sd(father))) %>%
  filter(z_father %in% -2:2) %>%
  ggplot() +  
  stat_qq(aes(sample = son)) +
  facet_wrap( ~ z_father)

Variance Explained

There is a link to the relevant section of the textbook: Variance explained

Key points

Conditioning on a random variable X can help to reduce variance of response variable Y.
The standard deviation of the conditional distribution is \(\mbox{SD}(Y \mid X=x) = \sigma_y\sqrt{1-\rho^2}\), which is smaller than the standard deviation without conditioning \(\sigma_y\).
Because variance is the standard deviation squared, the variance of the conditional distribution is \(\sigma_y^2(1-\rho^2)\).
In the statement “X explains such and such percent of the variability,” the percent value refers to the variance. The variance decreases by \(\rho^2\) percent.
The “variance explained” statement only makes sense when the data is approximated by a bivariate normal distribution.

There are Two Regression Lines

There is a link to the relevant section of the textbook: Warning: there are two regression lines

Key point

There are two different regression lines depending on whether we are taking the expectation of Y given X or taking the expectation of X given Y.

Code

# compute a regression line to predict the son's height from the father's height
mu_x <- mean(galton_heights$father)
mu_y <- mean(galton_heights$son)
s_x <- sd(galton_heights$father)
s_y <- sd(galton_heights$son)
r <- cor(galton_heights$father, galton_heights$son)
m_1 <-  r * s_y / s_x
b_1 <- mu_y - m_1*mu_x
m_1 # slope 1

## [1] 0.5027904

b_1 # intercept 1

## [1] 35.71249

# compute a regression line to predict the father's height from the son's height
m_2 <-  r * s_x / s_y
b_2 <- mu_x - m_2*mu_y
m_2 # slope 2

## [1] 0.4986676

b_2 # intercept 2

## [1] 33.96539

Assessment - Stratification and Variance Explained, Part 1

Look at the figure below. The slope of the regression line in this figure is equal to what, in words?

Scatter plot and regression line of son and father heights

[X] A. Slope = (correlation coefficient of son and father heights) * (standard deviation of sons’ heights / standard deviation of fathers’ heights)
[ ] B. Slope = (correlation coefficient of son and father heights) * (standard deviation of fathers’ heights / standard deviation of sons’ heights)
[ ] C. Slope = (correlation coefficient of son and father heights) / (standard deviation of sons’ heights * standard deviation of fathers’ heights)
[ ] D. Slope = (mean height of fathers) - (correlation coefficient of son and father heights * mean height of sons).

Why does the regression line simplify to a line with intercept zero and slope \(\rho\) when we standardize our x and y variables? Try the simplification on your own first!

[ ] A. When we standardize variables, both x and y will have a mean of one and a standard deviation of zero. When you substitute this into the formula for the regression line, the terms cancel out until we have the following equation: \(y_i = \rho x_i\).
[X] B. When we standardize variables, both x and y will have a mean of zero and a standard deviation of one. When you substitute this into the formula for the regression line, the terms cancel out until we have the following equation: \(y_i = \rho x_i\).
[ ] C. When we standardize variables, both x and y will have a mean of zero and a standard deviation of one. When you substitute this into the formula for the regression line, the terms cancel out until we have the following equation: \(y_i = \rho + x_i\).

What is a limitation of calculating conditional means?

[X] A. Each stratum we condition on (e.g., a specific father’s height) may not have many data points.
[X] B. Because there are limited data points for each stratum, our average values have large standard errors.
[X] C. Conditional means are less stable than a regression line.
[ ] D. Conditional means are a useful theoretical tool but cannot be calculated.

A regression line is the best prediction of Y given we know the value of X when:

[X] A. X and Y follow a bivariate normal distribution.
[ ] B. Both X and Y are normally distributed.
[ ] C. Both X and Y have been standardized.
[ ] D. There are at least 25 X-Y pairs.

Which one of the following scatterplots depicts an x and y distribution that is NOT well-approximated by the bivariate normal distribution?

[X] A.

[ ] B.

[ ] C.

[ ] D.

We previously calculated that the correlation coefficient between fathers’ and sons’ heights is 0.5.

Given this, what percent of the variation in sons’ heights is explained by fathers’ heights?

[ ] A. 0%
[X] B. 25%
[ ] C. 50%
[ ] D. 75%

When two variables follow a bivariate normal distribution, the variation explained can be calculated as \(\rho^2 \times 100\).

Suppose the correlation between father and son’s height is 0.5, the standard deviation of fathers’ heights is 2 inches, and the standard deviation of sons’ heights is 3 inches.

Given a one inch increase in a father’s height, what is the predicted change in the son’s height?

[ ] A. 0.333
[ ] B. 0.5
[ ] C. 0.667
[X] D. 0.75
[ ] E. 1
[ ] F. 1.5

The slope of the regression line is calculated by multiplying the correlation coefficient by the ratio of the standard deviation of son heights and standard deviation of father heights: \(\sigma_{son}/\sigma_{father}\).

Assessment - Stratification and Variance Explained, Part 2

In the second part of this assessment, you’ll analyze a set of mother and daughter heights, also from GaltonFamilies.

Define female_heights, a set of mother and daughter heights sampled from GaltonFamilies, as follows:

set.seed(1989, sample.kind="Rounding") #if you are using R 3.6 or later

## Warning in set.seed(1989, sample.kind = "Rounding"): non-uniform 'Rounding'
## sampler used

female_heights <- GaltonFamilies%>%     
    filter(gender == "female") %>%     
    group_by(family) %>%     
    sample_n(1) %>%     
    ungroup() %>%     
    select(mother, childHeight) %>%     
    rename(daughter = childHeight)

Calculate the mean and standard deviation of mothers’ heights, the mean and standard deviation of daughters’ heights, and the correlaton coefficient between mother and daughter heights.

Mean of mothers’ heights

mean(female_heights$mother)

## [1] 64.125

Standard deviation of mothers’ heights

sd(female_heights$mother)

## [1] 2.289292

Mean of daughters’ heights

mean(female_heights$daughter)

## [1] 64.28011

Standard deviation of daughters’ heights

sd(female_heights$daughter)

## [1] 2.39416

Correlation coefficient

cor(female_heights$mother, female_heights$daughter)

## [1] 0.3245199

Calculate the slope and intercept of the regression line predicting daughters’ heights given mothers’ heights. Given an increase in mother’s height by 1 inch, how many inches is the daughter’s height expected to change?

Slope of regression line predicting daughters’ height from mothers’ heights

r <- cor(female_heights$mother, female_heights$daughter)
s_y <- sd(female_heights$daughter)
s_x <- sd(female_heights$mother)
r * s_y/s_x

## [1] 0.3393856

Intercept of regression line predicting daughters’ height from mothers’ heights

mu_y <- mean(female_heights$daughter)
mu_x <- mean(female_heights$mother)
mu_y - (r * s_y/s_x)*mu_x

## [1] 42.51701

Change in daughter’s height in inches given a 1 inch increase in the mother’s height

r * s_y/s_x

## [1] 0.3393856

What percent of the variability in daughter heights is explained by the mother’s height?

r^2*100

## [1] 10.53132

A mother has a height of 60 inches.

What is the conditional expected value of her daughter’s height given the mother’s height?

m = r * s_y/s_x
b = mu_y - (r * s_y/s_x)*mu_x
x = 60
m*x+b

## [1] 62.88015

Section 2 - Linear Models Overview

In the Linear Models section, you will learn how to do linear regression.

After completing this section, you will be able to:

Use multivariate regression to adjust for confounders.
Write linear models to describe the relationship between two or more variables.
Calculate the least squares estimates for a regression model using the lm function.
Understand the differences between tibbles and data frames.
Use the do() function to bridge R functions and the tidyverse.
Use the tidy(), glance(), and augment() functions from the broom package.
Apply linear regression to measurement error models.

This section has four parts: Introduction to Linear Models, Least Squares Estimates, Tibbles, do, and broom, and Regression and Baseball.

Confounding: Are BBs More Predictive?

The textbook for this section is available here

Key points

Association is not causation!
Although it may appear that BB cause runs, it is actually the HR that cause most of these runs. We say that BB are confounded with HR.
Regression can help us account for confounding.

Code

# find regression line for predicting runs from BBs
bb_slope <- Teams %>% 
  filter(yearID %in% 1961:2001 ) %>% 
  mutate(BB_per_game = BB/G, R_per_game = R/G) %>% 
  lm(R_per_game ~ BB_per_game, data = .) %>% 
  .$coef %>%
  .[2]
bb_slope

## BB_per_game 
##   0.7353288

# compute regression line for predicting runs from singles
singles_slope <- Teams %>% 
  filter(yearID %in% 1961:2001 ) %>%
  mutate(Singles_per_game = (H-HR-X2B-X3B)/G, R_per_game = R/G) %>%
  lm(R_per_game ~ Singles_per_game, data = .) %>%
  .$coef  %>%
  .[2]
singles_slope

## Singles_per_game 
##        0.4494253

# calculate correlation between HR, BB and singles
Teams %>% 
  filter(yearID %in% 1961:2001 ) %>% 
  mutate(Singles = (H-HR-X2B-X3B)/G, BB = BB/G, HR = HR/G) %>%  
  summarize(cor(BB, HR), cor(Singles, HR), cor(BB,Singles))

##   cor(BB, HR) cor(Singles, HR) cor(BB, Singles)
## 1   0.4039313       -0.1737435      -0.05603822

Stratification and Multivariate Regression

The textbook for this section is available here

Key points

A first approach to check confounding is to keep HRs fixed at a certain value and then examine the relationship between BB and runs.
The slopes of BB after stratifying on HR are reduced, but they are not 0, which indicates that BB are helpful for producing runs, just not as much as previously thought.

Code

# stratify HR per game to nearest 10, filter out strata with few points
dat <- Teams %>% filter(yearID %in% 1961:2001) %>%
  mutate(HR_strata = round(HR/G, 1), 
         BB_per_game = BB / G,
         R_per_game = R / G) %>%
  filter(HR_strata >= 0.4 & HR_strata <=1.2)
  
# scatterplot for each HR stratum
dat %>% 
  ggplot(aes(BB_per_game, R_per_game)) +  
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm") +
  facet_wrap( ~ HR_strata)

## `geom_smooth()` using formula 'y ~ x'

# calculate slope of regression line after stratifying by HR
dat %>%  
  group_by(HR_strata) %>%
  summarize(slope = cor(BB_per_game, R_per_game)*sd(R_per_game)/sd(BB_per_game))

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 9 x 2
##   HR_strata slope
##       <dbl> <dbl>
## 1       0.4 0.734
## 2       0.5 0.566
## 3       0.6 0.412
## 4       0.7 0.285
## 5       0.8 0.365
## 6       0.9 0.261
## 7       1   0.511
## 8       1.1 0.454
## 9       1.2 0.440

# stratify by BB
dat <- Teams %>% filter(yearID %in% 1961:2001) %>%
  mutate(BB_strata = round(BB/G, 1), 
         HR_per_game = HR / G,
         R_per_game = R / G) %>%
  filter(BB_strata >= 2.8 & BB_strata <=3.9) 

# scatterplot for each BB stratum
dat %>% ggplot(aes(HR_per_game, R_per_game)) +  
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm") +
  facet_wrap( ~ BB_strata)

## `geom_smooth()` using formula 'y ~ x'

# slope of regression line after stratifying by BB
dat %>%  
  group_by(BB_strata) %>%
  summarize(slope = cor(HR_per_game, R_per_game)*sd(R_per_game)/sd(HR_per_game))

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 12 x 2
##    BB_strata slope
##        <dbl> <dbl>
##  1       2.8  1.52
##  2       2.9  1.57
##  3       3    1.52
##  4       3.1  1.49
##  5       3.2  1.58
##  6       3.3  1.56
##  7       3.4  1.48
##  8       3.5  1.63
##  9       3.6  1.83
## 10       3.7  1.45
## 11       3.8  1.70
## 12       3.9  1.30

Linear Models

The textbook for this section is available here

Key points

“Linear” here does not refer to lines, but rather to the fact that the conditional expectation is a linear combination of known quantities.
In Galton’s model, we assume \(Y\) (son’s height) is a linear combination of a constant and \(X\) (father’s height) plus random noise. We further assume that \(\epsilon_i\) are independent from each other, have expected value 0 and the standard deviation \(\sigma\) which does not depend on i.
Note that if we further assume that \(\epsilon\) is normally distributed, then the model is exactly the same one we derived earlier by assuming bivariate normal data.
We can subtract the mean from \(X\) to make \(\beta_0\) more interpretable.

Assessment: Introduction to Linear Models

When we stratified our regression lines for runs per game vs. bases on balls by the number of home runs, what happened?

[X] A. The slope of runs per game vs. bases on balls within each stratum was reduced because we removed confounding by home runs.
[ ] B. The slope of runs per game vs. bases on balls within each stratum was reduced because there were fewer data points.
[ ] C. The slope of runs per game vs. bases on balls within each stratum increased after we removed confounding by home runs.
[ ] D. The slope of runs per game vs. bases on balls within each stratum stayed about the same as the original slope.

We run a linear model for sons’ heights vs. fathers’ heights using the Galton height data, and get the following results:

> lm(son ~ father, data = galton_heights)

Call:
lm(formula = son ~ father, data = galton_heights)

Coefficients:
(Intercept)    father  
    35.71       0.50

Interpret the numeric coefficient for “father.”

[ ] A. For every inch we increase the son’s height, the predicted father’s height increases by 0.5 inches.
[X] B. For every inch we increase the father’s height, the predicted son’s height grows by 0.5 inches.
[ ] C. For every inch we increase the father’s height, the predicted son’s height is 0.5 times greater.

We want the intercept term for our model to be more interpretable, so we run the same model as before but now we subtract the mean of fathers’ heights from each individual father’s height to create a new variable centered at zero.

galton_heights <- galton_heights %>%
    mutate(father_centered=father - mean(father))

We run a linear model using this centered fathers’ height variable.

> lm(son ~ father_centered, data = galton_heights)

Call:
lm(formula = son ~ father_centered, data = galton_heights)

Coefficients:
(Intercept)    father_centered  
    70.45          0.50

Interpret the numeric coefficient for the intercept.

[X] A. The height of a son of a father of average height is 70.45 inches.
[ ] B. The height of a son when a father’s height is zero is 70.45 inches.
[ ] C. The height of an average father is 70.45 inches.

Suppose we fit a multivariate regression model for expected runs based on BB and HR:

\(E[R|BB=x_1, HR=x_2] = \beta_0+\beta_1x_1+\beta_2x_2\)

Suppose we fix \(BB = x_1\). Then we observe a linear relationship between runs and HR with intercept of:

[ ] A. \(\beta_0\)
[ ] B. \(\beta_0 + \beta_2x_2\)
[X] C. \(\beta_0 + \beta_1x_1\)
[ ] D. \(\beta_0 + \beta_2x_1\)

Which of the following are assumptions for the errors \(\epsilon_i\) in a linear regression model? Check ALL correct answers.

[X] A. The \(\epsilon_i\) are independent of each other
[X] B. The \(\epsilon_i\) have expected value 0
[X] C. The variance of \(\epsilon_i\) is a constant

Least Squares Estimates (LSE)

The textbook for this section is available here

Key points

For regression, we aim to find the coefficient values that minimize the distance of the fitted model to the data.
Residual sum of squares (RSS) measures the distance between the true value and the predicted value given by the regression line. The values that minimize the RSS are called the least squares estimates (LSE).
We can use partial derivatives to get the values for \(\beta_0\) and \(\beta_1\) in Galton’s data.

Code

# compute RSS for any pair of beta0 and beta1 in Galton's data
set.seed(1983)
galton_heights <- GaltonFamilies %>%
  filter(gender == "male") %>%
  group_by(family) %>%
  sample_n(1) %>%
  ungroup() %>%
  select(father, childHeight) %>%
  rename(son = childHeight)
rss <- function(beta0, beta1, data){
    resid <- galton_heights$son - (beta0+beta1*galton_heights$father)
    return(sum(resid^2))
}

# plot RSS as a function of beta1 when beta0=25
beta1 = seq(0, 1, len=nrow(galton_heights))
results <- data.frame(beta1 = beta1,
                      rss = sapply(beta1, rss, beta0 = 25))
results %>% ggplot(aes(beta1, rss)) + geom_line() + 
  geom_line(aes(beta1, rss))

The lm Function

The textbook for this section is available here

Key points

When calling the lm() function, the variable that we want to predict is put to the left of the ~ symbol, and the variables that we use to predict is put to the right of the ~ symbol. The intercept is added automatically.
LSEs are random variables.

Code

# fit regression line to predict son's height from father's height
fit <- lm(son ~ father, data = galton_heights)
fit

## 
## Call:
## lm(formula = son ~ father, data = galton_heights)
## 
## Coefficients:
## (Intercept)       father  
##     38.7646       0.4411

# summary statistics
summary(fit)

## 
## Call:
## lm(formula = son ~ father, data = galton_heights)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.4228 -1.7022  0.0333  1.5670  9.3567 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 38.76457    5.41093   7.164 2.03e-11 ***
## father       0.44112    0.07825   5.637 6.72e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.659 on 177 degrees of freedom
## Multiple R-squared:  0.1522, Adjusted R-squared:  0.1474 
## F-statistic: 31.78 on 1 and 177 DF,  p-value: 6.719e-08

LSE are Random Variables

The textbook for this section is available here

Key points

Because they are derived from the samples, LSE are random variables.
\(\beta_0\) and \(\beta_1\) appear to be normally distributed because the central limit theorem plays a role.
The t-statistic depends on the assumption that \(\epsilon\) follows a normal distribution.

Code

# Monte Carlo simulation
B <- 1000
N <- 50
lse <- replicate(B, {
  sample_n(galton_heights, N, replace = TRUE) %>% 
    lm(son ~ father, data = .) %>% 
    .$coef 
})
lse <- data.frame(beta_0 = lse[1,], beta_1 = lse[2,]) 

# Plot the distribution of beta_0 and beta_1
if(!require(gridExtra)) install.packages("gridExtra")

## Loading required package: gridExtra

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

library(gridExtra)
p1 <- lse %>% ggplot(aes(beta_0)) + geom_histogram(binwidth = 5, color = "black") 
p2 <- lse %>% ggplot(aes(beta_1)) + geom_histogram(binwidth = 0.1, color = "black") 
grid.arrange(p1, p2, ncol = 2)

# summary statistics
sample_n(galton_heights, N, replace = TRUE) %>% 
  lm(son ~ father, data = .) %>% 
  summary %>%
  .$coef

##               Estimate Std. Error  t value     Pr(>|t|)
## (Intercept) 34.4729422  8.6021831 4.007464 0.0002129225
## father       0.4990193  0.1240572 4.022493 0.0002030210

lse %>% summarize(se_0 = sd(beta_0), se_1 = sd(beta_1))

##       se_0      se_1
## 1 9.683973 0.1411404

Advanced Note on LSE

Although interpretation is not straight-forward, it is also useful to know that the LSE can be strongly correlated, which can be seen using this code:

lse %>% summarize(cor(beta_0, beta_1))

##   cor(beta_0, beta_1)
## 1          -0.9993386

However, the correlation depends on how the predictors are defined or transformed.

Here we standardize the father heights, which changes \(x_i\) to \(x_i - \bar{x}\).

B <- 1000
N <- 50
lse <- replicate(B, {
      sample_n(galton_heights, N, replace = TRUE) %>%
      mutate(father = father - mean(father)) %>%
      lm(son ~ father, data = .) %>% .$coef 
})

Observe what happens to the correlation in this case:

cor(lse[1,], lse[2,])

## [1] 0.1100929

Predicted Variables are Random Variables

The textbook for this section is available here

Key points

The predicted value is often denoted as \(\hat{Y}\), which is a random variable. Mathematical theory tells us what the standard error of the predicted value is.
The predict() function in R can give us predictions directly.

Code

# plot predictions and confidence intervals
galton_heights %>% ggplot(aes(son, father)) +
  geom_point() +
  geom_smooth(method = "lm")

## `geom_smooth()` using formula 'y ~ x'

# predict Y directly
fit <- galton_heights %>% lm(son ~ father, data = .) 
Y_hat <- predict(fit, se.fit = TRUE)
names(Y_hat)

## [1] "fit"            "se.fit"         "df"             "residual.scale"

# plot best fit line
galton_heights %>%
  mutate(Y_hat = predict(lm(son ~ father, data=.))) %>%
  ggplot(aes(father, Y_hat))+
  geom_line()

Assessment - Least Squares Estimates (LSE), part 1

The following code was used in the video to plot RSS with \(\beta_0=25\).

beta1 = seq(0, 1, len=nrow(galton_heights))
results <- data.frame(beta1 = beta1,
                      rss = sapply(beta1, rss, beta0 = 25))
results %>% ggplot(aes(beta1, rss)) + geom_line() + 
  geom_line(aes(beta1, rss), col=2)

In a model for sons’ heights vs fathers’ heights, what is the least squares estimate (LSE) for \(\beta_1\) if we assume \(\hat{\beta}_{0}\) is 36? Hint: modify the code above to do your analysis.

# compute RSS for any pair of beta0 and beta1 in Galton's data
set.seed(1983)
galton_heights <- GaltonFamilies %>%
  filter(gender == "male") %>%
  group_by(family) %>%
  sample_n(1) %>%
  ungroup() %>%
  select(father, childHeight) %>%
  rename(son = childHeight)
rss <- function(beta0, beta1, data){
    resid <- galton_heights$son - (beta0+beta1*galton_heights$father)
    return(sum(resid^2))
}

# plot RSS as a function of beta1 when beta0=36
beta1 = seq(0, 1, len=nrow(galton_heights))
results <- data.frame(beta1 = beta1,
                      rss = sapply(beta1, rss, beta0 = 36))
results %>% ggplot(aes(beta1, rss)) + geom_line() + 
  geom_line(aes(beta1, rss))

[ ] A. 0.65
[X] B. 0.5
[ ] C. 0.2
[ ] D. 12

The least squares estimates for the parameters \(\beta_0, \beta_1, \dots, \beta_n\) minimize the residual sum of squares.
Load the Lahman library and filter the Teams data frame to the years 1961-2001.

Run a linear model in R predicting the number of runs per game based on the number of bases on balls and the number of home runs. Remember to first limit your data to 1961-2001.

What is the coefficient for bases on balls?

fit <- Teams %>% filter(yearID %in% 1961:2001) %>% 
mutate(R_per_game = R / G, BB_per_game = BB / G, HR_per_game = HR / G) %>% 
lm(R_per_game ~ BB_per_game + HR_per_game, data = .)

summary(fit)

## 
## Call:
## lm(formula = R_per_game ~ BB_per_game + HR_per_game, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.87325 -0.24507 -0.01449  0.23866  1.24218 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.74430    0.08236   21.18   <2e-16 ***
## BB_per_game  0.38742    0.02701   14.34   <2e-16 ***
## HR_per_game  1.56117    0.04896   31.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3484 on 1023 degrees of freedom
## Multiple R-squared:  0.6503, Adjusted R-squared:  0.6496 
## F-statistic: 951.2 on 2 and 1023 DF,  p-value: < 2.2e-16

[X] A. 0.39
[ ] B. 1.56
[ ] C. 1.74
[ ] D. 0.027

The coefficient for bases on balls is 0.39; the coefficient for home runs is 1.56; the intercept is 1.74; the standard error for the BB coefficient is 0.027.

We run a Monte Carlo simulation where we repeatedly take samples of N = 100 from the Galton heights data and compute the regression slope coefficients for each sample:

B <- 1000
N <- 100
lse <- replicate(B, {
    sample_n(galton_heights, N, replace = TRUE) %>% 
    lm(son ~ father, data = .) %>% .$coef 
})

lse <- data.frame(beta_0 = lse[1,], beta_1 = lse[2,])

What does the central limit theorem tell us about the variables beta_0 and beta_1?

[X] A. They are approximately normally distributed.
[X] B. The expected value of each is the true value of \(\beta_0\) and \(\beta_1\) (assuming the Galton heights data is a complete population).
[ ] C. The central limit theorem does not apply in this situation.
[ ] D. It allows us to test the hypothesis that \(\beta_0 = 0\) and \(\beta_0 = 1\)

In an earlier video, we ran the following linear model and looked at a summary of the results.

$\beta_0$
> mod <- lm(son ~ father, data = galton_heights)
> summary(mod)

Call:
lm(formula = son ~ father, data = galton_heights)

Residuals:
   Min     1Q  Median     3Q    Max 
-5.902  -1.405  0.092    1.342  8.092 

Coefficients:
                 Estimate  Std. Error  t value     Pr(>|t|)  
(Intercept)     35.7125     4.5174       7.91    2.8e-13 ***
father           0.5028     0.0653       7.70    9.5e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
$\beta_0$

What null hypothesis is the second p-value (the one in the father row) testing?

[ ] A. \(\beta_1 = 1\), where \(\beta_1\) is the coefficient for the variable “father.”
[ ] B. \(\beta_1 = 0.503\), where \(\beta_1\) is the coefficient for the variable “father.”
[X] C. \(\beta_1 = 0\), where \(\beta_1\) is the coefficient for the variable “father.”

Which R code(s) below would properly plot the predictions and confidence intervals for our linear model of sons’ heights?

Select ALL that apply.

[ ] A.

          galton_heights %>% ggplot(aes(father, son)) +
    geom_point() +
    geom_smooth()

[X] B.

          galton_heights %>% ggplot(aes(father, son)) +
    geom_point() +
    geom_smooth(method = "lm")

[X] C.

          model <- lm(son ~ father, data = galton_heights)
predictions <- predict(model, interval = c("confidence"), level = 0.95)
data <- as.tibble(predictions) %>% bind_cols(father = galton_heights$father)

ggplot(data, aes(x = father, y = fit)) +
    geom_line(color = "blue", size = 1) + 
    geom_ribbon(aes(ymin=lwr, ymax=upr), alpha=0.2) + 
    geom_point(data = galton_heights, aes(x = father, y = son))

[ ] D.

          model <- lm(son ~ father, data = galton_heights)
predictions <- predict(model)
data <- as.tibble(predictions) %>% bind_cols(father = galton_heights$father)

ggplot(data, aes(x = father, y = fit)) +
    geom_line(color = "blue", size = 1) + 
    geom_point(data = galton_heights, aes(x = father, y = son))

Assessment - Least Squares Estimates, part 2

In Questions 7 and 8, you’ll look again at female heights from GaltonFamilies.

Define female_heights, a set of mother and daughter heights sampled from GaltonFamilies, as follows:

set.seed(1989, sample.kind="Rounding") #if you are using R 3.6 or later

## Warning in set.seed(1989, sample.kind = "Rounding"): non-uniform 'Rounding'
## sampler used

options(digits = 3)    # report 3 significant digits

female_heights <- GaltonFamilies %>%     
    filter(gender == "female") %>%     
    group_by(family) %>%     
    sample_n(1) %>%     
    ungroup() %>%     
    select(mother, childHeight) %>%     
    rename(daughter = childHeight)

Fit a linear regression model predicting the mothers’ heights using daughters’ heights.

What is the slope of the model?

fit <- lm(mother ~ daughter, data = female_heights)
fit$coef[2]

## daughter 
##     0.31

What the intercept of the model?

fit$coef[1]

## (Intercept) 
##        44.2

Predict mothers’ heights using the model.

What is the predicted height of the first mother in the dataset?

predict(fit)[1]

##    1 
## 65.6

What is the actual height of the first mother in the dataset?

female_heights$mother[1]

## [1] 67

We have shown how BB and singles have similar predictive power for scoring runs. Another way to compare the usefulness of these baseball metrics is by assessing how stable they are across the years. Because we have to pick players based on their previous performances, we will prefer metrics that are more stable. In these exercises, we will compare the stability of singles and BBs.

Before we get started, we want to generate two tables: one for 2002 and another for the average of 1999-2001 seasons. We want to define per plate appearance statistics, keeping only players with more than 100 plate appearances. Here is how we create the 2002 table:

bat_02 <- Batting %>% filter(yearID == 2002) %>%
    mutate(pa = AB + BB, singles = (H - X2B - X3B - HR)/pa, bb = BB/pa) %>%
    filter(pa >= 100) %>%
    select(playerID, singles, bb)

Now compute a similar table but with rates computed over 1999-2001. Keep only rows from 1999-2001 where players have 100 or more plate appearances, calculate each player’s single rate and BB rate per season, then calculate the average single rate (mean_singles) and average BB rate (mean_bb) per player over those three seasons.

How many players had a single rate mean_singles of greater than 0.2 per plate appearance over 1999-2001?

bat_99_01 <- Batting %>% filter(yearID %in% 1999:2001) %>%
    mutate(pa = AB + BB, singles = (H - X2B - X3B - HR)/pa, bb = BB/pa) %>%
    filter(pa >= 100) %>%
    group_by(playerID) %>%
    summarize(mean_singles = mean(singles), mean_bb = mean(bb))

## `summarise()` ungrouping output (override with `.groups` argument)

sum(bat_99_01$mean_singles > 0.2)

## [1] 46

How many players had a BB rate mean_bb of greater than 0.2 per plate appearance over 1999-2001?

sum(bat_99_01$mean_bb > 0.2)

## [1] 3

Use inner_join() to combine the `bat_02 table with the table of 1999-2001 rate averages you created in the previous question.

What is the correlation between 2002 singles rates and 1999-2001 average singles rates?

dat <- inner_join(bat_02, bat_99_01)

## Joining, by = "playerID"

cor(dat$singles, dat$mean_singles)

## [1] 0.551

What is the correlation between 2002 BB rates and 1999-2001 average BB rates?

cor(dat$bb, dat$mean_bb)

## [1] 0.717

Make scatterplots of mean_singles versus singles and mean_bb versus bb.

Are either of these distributions bivariate normal?

dat %>%
    ggplot(aes(singles, mean_singles)) +
    geom_point()

dat %>%
    ggplot(aes(bb, mean_bb)) +
    geom_point()

[ ] A. Neither distribution is bivariate normal.
[ ] B. singles and mean_singles are bivariate normal, but bb and mean_bb are not.
[ ] C. bb and mean_bb are bivariate normal, but singles and mean_singles are not.
[X] D. Both distributions are bivariate normal.

Fit a linear model to predict 2002 singles given 1999-2001 mean_singles.

What is the coefficient of mean_singles, the slope of the fit?

fit_singles <- lm(singles ~ mean_singles, data = dat)
fit_singles$coef[2]

## mean_singles 
##        0.588

Fit a linear model to predict 2002 bb given 1999-2001 mean_bb.

What is the coefficient of mean_bb, the slope of the fit?

fit_bb <- lm(bb ~ mean_bb, data = dat)
fit_bb$coef[2]

## mean_bb 
##   0.829

Assessment 8 - Advanced dplyr: Tibbles

What problem do we encounter when we try to run a linear model on our baseball data, grouping by home runs?

[ ] A. There is not enough data in some levels to run the model.
[X] B. The lm function does not know how to handle grouped tibbles.
[ ] C. The results of the lm function cannot be put into a tidy format.

Tibbles are similar to what other class in R?

[ ] A. Vectors
[ ] B. Matrices
[X] C. Data frames
[ ] D. Lists

Assessment 9 - Tibbles: Differences from Data Frames

What are some advantages of tibbles compared to data frames?

[X] A. Tibbles display better.
[X] B. If you subset a tibble, you always get back a tibble.
[X] C. Tibbles can have complex entries.
[X] D. Tibbles can be grouped.

Assessment 10 - do

What are two advantages of the do command, when applied to the tidyverse?

[ ] A. It is faster than normal functions.
[ ] B. It returns useful error messages.
[X] C. It understands grouped tibbles.
[X] D. It always returns a data.frame.

You want to take the tibble dat, which we’ve been using in this video, and run the linear model R ~ BB for each strata of HR. Then you want to add three new columns to your grouped tibble: the coefficient, standard error, and p-value for the BB term in the model.

You’ve already written the function get_slope, shown below.

get_slope <- function(data) {
  fit <- lm(R ~ BB, data = data)
  sum.fit <- summary(fit)

  data.frame(slope = sum.fit$coefficients[2, "Estimate"], 
             se = sum.fit$coefficients[2, "Std. Error"],
             pvalue = sum.fit$coefficients[2, "Pr(>|t|)"])
}

What additional code could you write to accomplish your goal?

[ ] A.

dat %>% 
  group_by(HR) %>% 
  do(get_slope)

[X] B.

dat %>% 
  group_by(HR) %>% 
  do(get_slope(.))

[ ] C.

dat %>% 
  group_by(HR) %>% 
  do(slope = get_slope(.))

[ ] D.

dat %>% 
  do(get_slope(.))

Assessment 11 - broom

The output of a broom function is always what?

[X] A. A data.frame
[ ] B. A list
[ ] C. A vector

You want to know whether the relationship between home runs and runs per game varies by baseball league. You create the following dataset:

dat <- Teams %>% filter(yearID %in% 1961:2001) %>%
  mutate(HR = HR/G,
         R = R/G) %>%
  select(lgID, HR, BB, R)

What code would help you quickly answer this question?

[X] A.

dat %>% 
  group_by(lgID) %>% 
  do(tidy(lm(R ~ HR, data = .), conf.int = T)) %>% 
  filter(term == "HR")

[ ] B.

dat %>% 
  group_by(lgID) %>% 
  do(glance(lm(R ~ HR, data = .)))

[ ] C.

dat %>% 
  do(tidy(lm(R ~ HR, data = .), conf.int = T)) %>% 
  filter(term == "HR")

[ ] D.

dat %>% 
  group_by(lgID) %>% 
  do(mod = lm(R ~ HR, data = .))

Assessment 12 - Building a Better Offensive Metric for Baseball

What is the final linear model we use to predict runs scored per game?

[ ] A. lm(R ~ BB + HR)
[ ] B. lm(HR ~ BB + singles + doubles + triples)
[X] C. lm(R ~ BB + singles + doubles + triples + HR)
[ ] D. lm(R ~ singles + doubles + triples + HR)

We want to estimate runs per game scored by individual players, not just by teams. What summary metric do we calculate to help estimate this?

Look at the code from the video for a hint:

pa_per_game <- Batting %>% 
  filter(yearID == 2002) %>% 
  group_by(teamID) %>%
  summarize(pa_per_game = sum(AB+BB)/max(G)) %>% 
  .$pa_per_game %>% 
  mean

[ ] A. pa_per_game: the mean number of plate appearances per team per game for each team
[ ] B. pa_per_game: the mean number of plate appearances per game for each player
[X] C. pa_per_game: the number of plate appearances per team per game, averaged across all teams

Imagine you have two teams. Team A is comprised of batters who, on average, get two bases on balls, four singles, one double, and one home run. Team B is comprised of batters who, on average, get one base on balls, six singles, two doubles, and one triple.

Which team scores more runs, as predicted by our model?

[ ] A. Team A
[X] B. Team B
[ ] C. Tie
[ ] D. Impossible to know

Assessment 13 - On Base Plus Slugging (OPS)

The on-base-percentage plus slugging percentage (OPS) metric gives the most weight to:

[ ] A. Singles
[ ] B. Doubles
[ ] C. Triples
[X] D. Home Runs

Assessment 14 - Regression Fallacy

What statistical concept properly explains the “sophomore slump”?

[X] A. Regression to the mean
[ ] B. Law of averages
[ ] C. Normal distribution

Assessment 15 - Measurement Error Models

In our model of time vs. observed_distance, the randomness of our data was due to:

[ ] A. sampling
[ ] B. natural variability
[X] C. measurement error

Which of the following are important assumptions about the measurement errors in this experiment?

[X] A. The measurement error is random
[X] B. The measurement error is independent
[X] C. The measurement error has the same distribution for each time i

Which of the following scenarios would violate an assumption of our measurement error model?

[ ] A. The experiment was conducted on the moon.
[X] B. There was one position where it was particularly difficult to see the dropped ball.
[ ] C. The experiment was only repeated 10 times, not 100 times.

Section 3 - Confounding Overview

In the Confounding section, you will learn what is perhaps the most important lesson of statistics: that correlation is not causation.

After completing this section, you will be able to:

Identify examples of spurious correlation and explain how data dredging can lead to spurious correlation.
Explain how outliers can drive correlation and learn to adjust for outliers using Spearman correlation.
Explain how reversing cause and effect can lead to associations being confused with causation.
Understand how confounders can lead to the misinterpretation of associations.
Explain and give examples of Simpson’s Paradox.

This section has one part: Correlation is Not Causation.

The textbook for this section is available here

Assessment 1 - Correlation is Not Causation: Spurious Correlation

In the video, we ran one million tests of correlation for two random variables, X and Y.

How many of these correlations would you expect to have a significant p-value (p>0.05), just by chance?

[ ] A. 5,000
[X] B. 50,000
[ ] C. 100,000
[ ] D. It’s impossible to know

Which of the following are examples of p-hacking?

[X] A. Looking for associations between an outcome and several exposures and only reporting the one that is significant.
[X] B. Trying several different models and selecting the one that yields the smallest p-value.
[X] C. Repeating an experiment multiple times and only reporting the one with the smallest p-value.
[ ] D. Using a Monte Carlo simulations in an analysis.

Assessment 2 - Correlation is Not Causation: Outliers

The Spearman correlation coefficient is robust to outliers because:

[ ] A. It drops outliers before calculating correlation.
[ ] B. It is the correlation of standardized values.
[X] C. It calculates correlation between ranks, not values.

Assessment 3 - Correlation is Not Causation: Reversing Cause and Effect

Which of the following may be examples of reversed cause and effect?

[X] A. Past smokers who have quit smoking may be more likely to die from lung cancer.
[ ] B. Tall fathers are more likely to have tall sons.
[X] C. People with high blood pressure tend to have a healthier diet.
[X] D. Individuals in a low social status have a higher risk of schizophrenia.

Assessment 4 - Correlation is Not Causation: Confounders

What can you do to determine if you are misinterpreting results because of a confounder?

[ ] A. Nothing, if the p-value says the result is significant, then it is.
[X] B. More closely examine the results by stratifying and plotting the data.
[ ] C. Always assume that you are misinterpreting the results.
[ ] D. Use linear models to tease out a confounder.

Look again at the admissions data using ?admissions. What important characteristic of the table variables do you need to know to understand the calculations used in this video? Select the best answer.

[ ] A. The data is from 1973.
[ ] B. The columns “major” and “gender” are of class character, while “admitted” and “applicants” are numeric.
[ ] C. The data is from the “dslabs” package.
[X] D. The column “admitted” is the percent of student admitted, while the column “applicants” is the total number of applicants.

In the example in the video, major selectivity confounds the relationship between UC Berkley admission rates and gender because:

[ ] A. It was harder for women to be admitted to UC Berkeley.
[X] B. Major selectivity is associated with both admission rates and with gender, as women tended to apply to more selective majors.
[ ] C. Some majors are more selective than others
[ ] D. Major selectivity is not a confounder.

Assessment 5 - Simpson’s Paradox

Admission rates at UC Berkeley are an example of Simpson’s Paradox because:

[X] A. It appears that men have higher a higher admission rate than women, however, after we stratify by major, we see that on average women have a higher admission rate than men.
[ ] B. It was a paradox that women were being admitted at a lower rate than men.
[ ] C. The relationship between admissions and gender is confounded by major selectivity.

Data Science Linear Regression

Learning Objectives

Course Overview

Introduction to Linear Regression

Linear Models

Confounding

Section 1 - Introduction to Regression Overview

Motivating Example: Moneyball

Baseball basics

Bases on Balls or Stolen Bases?

Assessment - Baseball as a Motivating Example

Correlation

Correlation Coefficient

Sample Correlation is a Random Variable

Assessment - Correlation

Anscombe’s Quartet/Stratification

Bivariate Normal Distribution

Variance Explained

There are Two Regression Lines

Assessment - Stratification and Variance Explained, Part 1

Assessment - Stratification and Variance Explained, Part 2

Section 2 - Linear Models Overview

Confounding: Are BBs More Predictive?

Stratification and Multivariate Regression

Linear Models

Assessment: Introduction to Linear Models

Least Squares Estimates (LSE)

The lm Function

LSE are Random Variables

Advanced Note on LSE

Predicted Variables are Random Variables

Assessment - Least Squares Estimates (LSE), part 1

Assessment - Least Squares Estimates, part 2

Assessment 8 - Advanced dplyr: Tibbles

Assessment 9 - Tibbles: Differences from Data Frames

Assessment 10 - do

Assessment 11 - broom

Assessment 12 - Building a Better Offensive Metric for Baseball

Assessment 13 - On Base Plus Slugging (OPS)

Assessment 14 - Regression Fallacy

Assessment 15 - Measurement Error Models

Section 3 - Confounding Overview

Assessment 1 - Correlation is Not Causation: Spurious Correlation

Assessment 2 - Correlation is Not Causation: Outliers

Assessment 3 - Correlation is Not Causation: Reversing Cause and Effect

Assessment 4 - Correlation is Not Causation: Confounders

Assessment 5 - Simpson’s Paradox